#include <matrix.hpp>
#include <random.hpp>

#include <algorithm>
#include <iostream>

#include <iomanip>

using namespace std;
using namespace sm;

int main()
{
	cout.precision(1);

	dynamic_matrix<double> a( 10, 10 );
	feng::variate_generator<double> vg;
	copy( vg.begin(), vg.begin()+100, a.begin() );

	lu_decomposition<dynamic_matrix<double> > ld( a );
	dynamic_matrix<double> l = ld.lower();
	dynamic_matrix<double> u = ld.upper();
	dynamic_matrix<double> lu = l*u;
	cout 	<< "\na is:\n" << a
			<< "\nl is:\n" << l
			<< "\nu is:\n" << u 
			<< "\nl*u is:\n" << lu
			<< "\na-lu is:\n" << a-lu
			<< "\n";
	static_matrix<double, 10, 10> b;
	copy( vg.begin(), vg.begin()+100, b.begin() );

	lu_decomposition<static_matrix<double, 10, 10> > lsd(b);
	static_matrix<double, 10, 10> sl = lsd.lower();
	static_matrix<double, 10, 10> su = lsd.upper();
	static_matrix<double, 10, 10> slu = sl*su;
	cout 	<< "\nb is:\n" << b
			<< "\nsl is:\n" << sl
			<< "\nsu is:\n" << su 
			<< "\nsl*su is:\n" << slu
			<< "\nb-slu is:\n" << b-slu
			<< "\n";
	return 0;
}
